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Abstract 


Nonlinear entropy stability and a summation-by-parts framework are used to derive provably stable, 
polynomial-based spectral collocation methods of arbitrary order. The new methods are closely 
related to discontinuous Galerkin spectral collocation methods commonly known as DGFEM, but 
exhibit a more general entropy stability property. Although the new schemes are applicable to a 
broad class of linear and nonlinear conservation laws, emphasis herein is placed on the entropy 
stability of the compressible Navier-Stokes equations. 
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1 Introduction 

A current organizational research goal of NASA’s Revolutionary Computational AeroSciences sub- 
project of the Aeronautical Sciences Project is to develop next generation high-order numerical 
algorithms for use in large eddy simulations (LES) and hybrid Reynolds-averaged Navier- Stokes 
(RANS)-LES simulations of complex separated flow. These algorithms must be suitable for sim- 
ulations of highly nonlinear next generation turbulence models across the subsonic, transonic and 
supersonic speed regimes. 

Although high-order techniques are well suited for LES, most lack robustness when the solution 
contains discontinuities or even under-resolved physical features. Although a variety of mathemat- 
ically rigorous stabilization techniques have been developed for second-order methods (e.g., total 
variation diminishing (TVD) limiters [1], and entropy stability [2]), extending these techniques to 
high-order formulations has been problematic. A typical consequence is loss of design order accu- 
racy at local extrema or insufficient stabilization. It is possible by using essentially nonoscillatory 
(ENO) [3,4] and weighted ENO (WENO) [5,6] schemes, to achieve high-order design accuracy 
away from captured discontinuities, and maintain sharp “nearly monotone” captured shocks. Un- 
fortunately, nonoscillatory schemes experience instabilities in less than ideal circumstances (i.e. , 
curvilinear mapped grids or expansion of flows into vacuum). Because these schemes are largely 
based on stencil biasing heuristics rather than mathematical stability proofs and the theory that 
does exist is not sharp [7,8], there is little to guide further development efforts focused on alleviating 
the instabilities; that is, until recently. 

Fisher and Carpenter [9,10] provide a general procedure for developing entropy conservative and 
entropy stable schemes of any order, for broad classes of spatial operators. The work generalizes 
entropy stability results appearing in the literature by several authors over the past two decades. A 
brief overview of the evolutionary developments of entropy stability is now presented. Nearly three 
decades ago, entropy conservative schemes that discretely satisfy an entropy conservation property 
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are constructed by Tadmor for second-order finite volume methods [2, 11]. These schemes were 
extended to high-order periodic domains by LeFloch and Rhode [12]. Finding a computationally 
friendly discrete entropy flux was a major obstacle that was alleviated recently for the Navier-Stokes 
equations through the work of Ismail and Roe [13]. A methodology for constructing entropy stable 
schemes satisfying a cell entropy inequality and capable of simulating flows with shocks in periodic 
domains was developed by Fjordholm et al. [14]. Recently, Fisher and Carpenter [9,10] present 
multi-domain proofs of entropy conservation and stability based on diagonal norm summation-by- 
parts (SBP) operators, that yield entropy stable methods on finite domains. Generalization to 
arbitrary Cartesian domains follows immediately using simultaneous-approximation-term (SAT) 
penalty type interface conditions [15] between adjoining domains. 

Although the primary focus of references [9,10] is entropy stable WENO finite-difference schemes, 
all proofs immediately generalize to a broad class of summation-by-parts (SBP-SAT) operators. 
Indeed, any spatial discretization that may be expressed as a non-dissipative, diagonal norm, SBP- 
SAT operator can be implemented in an entropy conservative and stable fashion. Thus, it is now 
possible to construct entropy conservative and stable formulations of arbitrary order for many 
popular discrete operators. 

Spectral collocation operators are readily expressed in SBP-SAT form [16-19], although not all 
may be expressed as diagonal norm SBP operators; e.g., the mass matrix of a Chebyshev operator 
is full. Legendre collocation schemes, however, may be expressed as diagonal norm SBP operators, 
and therefore satisfy the sufficient conditions for an entropy stable implementation. Thus the focus 
herein is developing conservation form, entropy conservative, Legendre spectral collocation schemes 
for the Navier-Stokes equations. 1 

To place this work into the appropriate context, note that the SBP-SAT entropy stability proofs 
developed herein (and references [9,10]), are not the first appearing in the high-order literature, 
although they do have some unique properties. Nonlinear analysis for the Navier-Stokes equa- 
tions appears in the work of Hughes et al. [20] in the context of Galerkin and Petrov-Galerkin 
finite-element methods (FEM). Entropy stability of the FEM follows immediately by rotating the 
conservative equations into symmetric form followed by a conventional FEM implementation. Sym- 
metrizing the equations, however, raises the question of whether the method is consistent with the 
Lax-Wendroff theorem. The SBP-SAT stability proofs are implemented in conservation form and 
do not have this ambiguity. 

Likewise, entropy stability proofs for alternative nonlinear equations appear in many finite- 
element texts, e.g., see Hesthaven and Warburten [21] for a discussion of Burgers’ equation. Exten- 
sion to the compressible Navier-Stokes equations in conservation form, has not been forthcoming 
to our knowledge. Indeed, a fundamental obstacle in FEM proofs is the requirement for exact 
integration formulae, a feat that is all but impossible for the compressible Navier-Stokes equations. 
(Recasting the equations in entropy variables is not an option because of inconsistency with the 
Lax-Wendroff theorem.) Again, the SBP-SAT entropy stability proofs do not suffer these limita- 
tions. 

The entropy stable discrete operators developed herein are an important step towards a provably 
stable simulation methodology of arbitrary order for complex geometries. All proofs generalized 
immediately to 3D via tensor product arithmetic. The extension of the entropy stable methods 

1 An entropy stable correction for these Legendre spectral collocation scheme follows immediately from the results 
in [9,10] (e.g., a dissipative shock capturing numerical method such as WENO), but this development effort is deferred 
to a later paper. 


3 


to generalized (3D) curvilinear coordinates and multi-domain configurations is included in a com- 
panion paper. Two major hurdles remain on the path towards L 2 stability of the compressible 
Navier-Stokes equations. First, well-posed physical boundary conditions are needed that preserve 
the entropy stability property of the interior operator. Second, a fully discrete operator is needed 
that preserves the semi-discrete entropy stability properties (e.g., the trapezoidal rule) [11], and 
maintains positivity of the density and temperature. 

The organization of this paper is as follows. The theory of SBP-SAT operators and their rela- 
tionship to polynomial spectral collocation formulations, is presented in Section 2. This discussion 
is tutorial in extent, and may be skipped by readers familiar with SBP-SAT nomenclature and op- 
erators. Section 3 presents an introduction to continuous entropy analysis followed by semi-discrete 
analysis that demonstrates the entropy- mimetic properties of diagonal norm SBP operators. The 
analysis is valid for arbitrarily high-order accurate Legendre spectral collocation operators. Section 
4 presents inviscid and viscous coupling conditions used to connect adjoining elements. Details of 
the implementation of entropy conservative operators in the context of the compressible Euler and 
Navier-Stokes equations are provided in Section 5. Finally, the accuracy of the resulting high-order 
schemes are demonstrated in Section 6, and conclusions are discussed in Section 7. 

2 Methodology 

Consider the calorically perfect Navier-Stokes equations, which may be expressed in the form 

qt + (f)xi = xen, f<E[0,oo), 

Bq = gb, x G <9D, te[0,oo), (2.1) 

q(x,0) = go (x), ieH, 

where the Cartesian coordinates, x = (aq, aq, x^) T , and time, t. are independent variables, and 
index sums are implied. The vectors q, /*, and /bd* are the conserved variables, the conserved 
inviscid fluxes and viscous fluxes, respectively. Without loss of generality, a three dimensional box 

n = [®1 ,xf] x [x%,x%] X [a$,xf] 

is chosen as our computational domain with dVL representing the boundary of the domain. The 
boundary vector, , is assumed to contain well-posed Dirichlet /Neumann data. We have omit- 
ted a detailed description of the three-dimensional Navier-Stokes equations, which may be found 
elsewhere [22]. 

2.1 Summation-by-parts Operators 

2.1.1 First Derivative 

First derivative operators that satisfy the summation-by-parts (SBP) convention, discretely mimic 
the integration-by-parts condition 


XR 

J (pq x dx 

X L 



( 2 . 2 ) 
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This mimetic property is achieved by constructing the first derivative approximation, T>4>, with an 
operator in the form 

v = v- l Q, V = V T , C T VC>0, c/o, 

Q t = B-Q, B = diag (—1, 0, . . . , 0, 1) . 

While it is not true in general that V is diagonal, herein the focus is exclusively on diagonal norm 
SBP operators, based on fixed element-based polynomials. The matrix V may be thought of as a 
mass matrix in the context of Galerkin finite-elements (FEM), incorporates the local grid spacing 
into the derivative definition. The nearly skew-symmetric matrix, Q, is an undivided differencing 
operator where all rows sum to zero and the first and last column sum to —1 and 1, respectively. 
The accuracy of the first derivative operator, T>, may be expressed as 

4> x {-x) =V(f) + T( yP+ i), (2.4) 

where T(p+\) is the truncation error of the approximation, and p is the order of the polynomial. 
Integration in the approximation space is conducted using an inner product with the appropriate 
integration weights contained in the norm V , 

XR 

J 4>q x dx « 4> t W(\, cj) = (0(xi),^(.t 2 ), ■ ■ • ,4>{x n )) T . (2.5) 

XL 

Using the definition in 2.3, the SBP property is demonstrated, 

4> T VV^ l Qc{ = 4 > t (B- Q t ) q = 4> N qN - faqi - 4> T v T v q. ( 2 . 6 ) 

The specific operators used in this work is shown in Appendix A. 

2.1.2 The Second Derivative 

The viscous approximations are written in general as 

(^x)q x ^)) x =V 2 ^) q + r p ^\ (2.7) 

also satisfy the SBP condition. Integration by parts yields 

xr xr 

J cj)('&q x ) x dx = </y&q x \l* ~ J 4>x$qxdx. (2.8) 

XL X L 

The second derivative variable coefficient operator resulting from two applications of the first deriva- 
tive may be manipulated for diagonal norm, V, into the expression 

W) = (-V t T[V\D + B[ti\V) , V T V[d]V = (V t V[$}V) T , [0] = diag (0(x)) , 

C T (v T v[$]v) c > o, C T [tf]C>o, vc. 

The P-norm inner product yields the expression 

< t> T VV _1 (-V T V[V]V + B[ti\D) q = 4> T B[$\D q - 0 T {V T V[d]V) q (2.10) 
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which is the form used to show stability of the viscous terms. It is clear that the continuous interface 
terms are mimicked. Likewise, based on the definition 2.9 the expression, 



follows immediately. 


2.1.3 Complementary Grids 

Most existing entropy analysis is performed in indicial notation on a staggered set of solution and 
flux points. For example, Tadmor’s telescoping entropy flux relation (fully defined in section 3.2.2) 
is written as 


and relates solution point data Wi, Wi+i, ifi, V’i+i with a flux fl located between the grid points. 
Conventional SBP operators are not directly applicable to this form of analysis; generalized opera- 
tors suitable for a staggered grid implementation are now developed. The new operators satisfy a 
generalized summation-by-parts property. 

Define on the interval — 1 < x < 1, the vectors of discrete solution points 


control volume edges employed in the finite volume method. The distribution of the flux points 
depends on the discretization operator. The spacing between the flux points is implicitly contained 


(m+ 1 - Wi) T fi = 4>i + 1 - 4>i, 



( 2 . 11 ) 


Since the approximate solution is constructed at these points, they are denoted the solution points. 
It is useful to define a set of intermediate points prescribing bounding control volumes about each 
solution point. These (N + 1) points are denoted flux points as they are similar in nature to the 


in the norm V: the diagonal elements of V are equal to the spacing between flux points, 


x = (xo, xi, . . . xn) T , xo = xi, xn = xn, 
Xi - Xi-i = i = 1,2,..., IV. 


( 2 . 12 ) 


In operator notation, this is equivalent to 


Ax = VI ; 1 = (1,1,...,1) T 


(2.13) 


where A is as defined in equation 2.16. Note that in 2.12, the first and last flux points are coincident 
with the first and last solution points, which enables the endpoint fluxes to be consistent, 


/o = f(qi), /tv = / (qn) ■ 


(2-14) 


This duality is needed to define unique operators and is important in proving entropy stability. 
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2.1.4 Telescopic Flux Form 


All SBP derivative operators, T>, can be manipulated into the telescopic flux form, 

fM) = r~ 1 Qi + T {p+1) = V-'M + T^y 

where the N x (N + 1) matrix A is 


A = 


( -1 
0 

0 
0 
0 


1 

-1 

0 

0 

0 


0 

1 

0 

0 


0 0 0 \ 

0 0 0 

■••0 0 . 
-1 1 0 

0 -11/ 


(2.15) 


(2.16) 


that calculates the undivided difference of the two adjacent fluxes. All conservative and accurate 
flux gradients may be constructed in the form of 2.15 for all SBP operators, Q, a fact that is 
reiterated in the following lemma presented without proof. (The original proof appears in reference 
[23]). 

Lemma 2.1. All differentiation matrices that satisfy the SBP convention given in eq. (2.3) are 
telescoping operators in the norm V . 


This telescopic flux form admits a generalized SBP property. All SBP operators defined in 
equation 2.3 can be manipulated to transfer the action of the discrete derivative onto a test function 
with an equivalent order of approximation. The telescopic flux form defined in equation 2.15 
combined with the flux consistency condition results in a more generalized relation, 

ct> T VV ~ l Af = 4> t (B - A)f = f(q N )f N - /(gi)0i - 0Af, 


(2.17) 


where 



f o 

-1 

0 

0 

0 

0 ^ 


( 

-1 

0 

0 

0 

0 

0 \ 


0 

1 

-1 

0 

0 

0 



0 

0 

0 

0 

0 

0 

A = 

0 

0 



0 

0 

, B = 


0 

0 



0 

0 


0 

0 

0 

1 

-1 

0 



0 

0 

0 

0 

0 

0 


V o 

0 

0 

0 

1 

0 / 


V 

0 

0 

0 

0 

0 

1 / 


and 


^ T A = ^ + 0(N~ 1 ) 


This is equivalent to the commonly used explanation of summation-by-parts in indicial form, 

iV-l 

! 4>i ( fi - fi- 1) = f(QN)(f>N - f(qi)h -£/< *+! ~ &) • ( 2 - 18 ) 


N 

£■ 

i — 1 


i = 1 


The action of the derivative is still moved onto the test function but at first order accuracy. Note 
that although this generalized property is used herein to construct entropy conservative fluxes, it 
is also instrumental for satisfying the Lax-Wendroff theorem [24] in weak form. 

Likewise, the variable coefficient viscous operators presented in Section 2.1.2 may be expressed 
in the form 

OMx))* ~ (~V T V[d]V + B[V]V) q = (2.19) 

and satisfy a telescoping conservation property which is identical to that of the inviscid terms. 
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2.1.5 The Semi-Discrete Operator 


Based on the previous discussion of SBP operators and their equivalent telescoping form, the semi- 
discrete form of equation 2.1 becomes 

q* = + 'DMiPi q) + r- l gb , = v- l /\ r (-f + f (,;)i ) + v- 1 ®,, 

q(x,0) = g 0 (x), ilSl, 

with g b containing the enforcement of boundary conditions. Full implementation details, including 
the viscous Jacobian [c] - tensors are included in previous works [23,25] and elsewhere [26-30]. 

Remark. It is not necessary to implement an SBP scheme in flux form, but is the natural form 
to add dissipation while retaining consistency with the Lax Wendroff theorem [23]. Furthermore, 
the senri-discrete entropy analysis presented in Section 5.1.3 relies on the existence of the flux form. 


2.2 Spectral Discretization Operators 


Spectral collocation methods are commonly implemented on computational grids based on the nodes 
of Gauss-quadrature formulas (i.e. , Gauss, Gauss Radau, or Gauss Lobatto (GL)). These smooth 
but nonuniform grids are highly clustered at the boundaries of the domain, in stark contrast to the 
uniform grids used in conventional finite difference methods. 

The numerical methods developed herein are all collocated at the Legendre Gauss-Labotto 
Legendre (LGL) points, and include both end points of the interval. This distribution that includes 
the end points, allows the operators to be written in terms of flux differences, analogous to a finite 
volume method and consistent with equations 2.17 and 2.19. The complete discretization operator 
for the p = 4 element is illustrated in Figure 1. 


fo fi 
fi 
Y o 

j|c 

Xl 


h 

Ui 

X-2 


h 


i 


fz 

U 2 
2-' 3 


h 


1 


£o Xl 


X2 


^3 

-1 = 1 
1 10 

[3 

-16 

0 

+ 16 

V 7 

45 

45 


h 

U3 

Xi 


fi h 
/s 

U4 

>|( II 

X 5 
Xi x 5 


+9 + 1 
10 ' 


Figure 1. The one-dimensional discretization for p = 4 Legendre collocation is illustrated. Solution 
points are denoted by • and flux points are denoted by X . 


2.2.1 Lagrange Polynomials 

Define the Lagrange polynomials on the discrete points x as 


L\{x) 


(1 —x)L(x) 

ai(-i) 


L(x) = ( x 

i f , n{x) 


- X 2 )(x - X 3 )...(X - X n - 2 )(x - X N -i) 

_ . r .ci _ . 

2Z/(+l) ’ ^ (x—Xj)L'(xj) 
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2<j<N-l 


(2.21) 


With a slight abuse of notation, define the vector of Lagrange polynomials as 


L(.x) = [Li(x),L 2 (x), • • • ,L n _i(x),L n (x)] T 


( 2 . 22 ) 


2.3 Differentiation 


Assume that a smooth and (infinitely) differentiable function, f(x), is defined on the interval 
— 1 < x < 1. Reading the function / and derivative f at the discrete points, x, yields the vectors 


f( x ) = [f(xi),f(x 2 ),---,f{x N -i),f(x N )] T ; 

f/ ( x ) = 


The interpolation polynomial, /tv(.t), that collocates /(x) at the points x is given by the con- 
traction 

f(x) ~ f(N- 1)0*0 = [L(x)] T f( x )- ( 2 - 24 ) 

Derivative operators expressed in terms of the Lagrange polynomials on the interval are derived in 
the following theorem, presented without proof. (The proof appears in many texts, e.g., reference 
[31].) 

Theorem 2.2. The derivative operator that exactly differentiates an arbitrary n th order polynomial 
at the collocation points, x, is 

V = [L’jixif . (2.25) 

The elements ofV are dij for 1 <i,j<n . 

An equivalent representation of the differentiation operator may also be used, that satisfies 
all the requirements for being an SBP operator (but in general will not be a diagonal norm SBP 
operator). 

Theorem 2.3. The derivative operator that exactly differentiates an arbitrary p th order polynomial 
(p = N — 1) at the collocation points, x, may be expressed as 

V = V~ 1 Q. (2.26) 

Proof. First note that in addition to expression (2.33), the exact derivative c -^~ of the function 
f(x) may be approximated by 

f(x) « = [L(x)] T f'(x) . (2.27) 

The Galerkin statement demands that the integral error between the two expressions be orthogonal 
to the basis set which in this case are the Lagrange polynomials L(x). This statement may be 
expressed as 

i 

/ L(x) |^[L(x)] T f / (x) — [L (x)] f(x)^ dx = 0 , (2.28) 

-l 

or in the equivalent form 

Rf'(x) = Qf(x), (2.29) 


9 


with 


V = f* 1 L(x)[L(x)] 7 dx ; Q = L(x)[L/(x)] T dx . (2.30) 

Equation (2.26) follows immediately when V is symmetric positive definite (SPD), and therefore 
invertible. 

The symmetry of V follows immediately from definition (2.30). Positive definiteness of V is 
established by pre- and post-multiply V by an arbitrary nonzero discrete vector, -0, which yields 
the expression 

ip't'V'ip = J^_ 1 'iJj t L(x)[L(x)] 1 ip dx = f^_ 1 , tp(x) 2 dx (2-31) 

that is strictly greater than zero, unless ij} is the null vector. Thus, the matrix V is SPD, therefore 
invertible, and (2.26) follows immediately. □ 

Remark Once the invertibility of V established, equation (2.26) may be proven directly by 
showing that V V = Q. 

A proof that Q is nearly skew-symmetric is as follows. 

Theorem 2.4. The matrix Q = /j L 1 L(.t)[L ' (x)] T dx is structurally of the form 

Q + Q t = B. (2.32) 

Thus, by virtue of the structure of V and Q, the differentiation operator, V, is indeed an SBP 
operator defined by 2.3. 

Proof. Integrating by parts the definition of Q yields the expression 

l l 

Q = J L(x)[L (x)] dx = L(x)(+1)[L(x)(+ 1)] T — L(.x)(— l)[L(x)(-l)] T - J L (x)[L(x)] T dx. 

-l -l 

(2.33) 

All Lagrange polynomials based on the Gauss-Labotto collocation points vanish on the boundaries 
for 1 < i,j < N. Thus, the boundary matrices reduce to the form 

L(x)(+1)[L(x)(+1)] t - L(x)(-l)[L(x)(-l)] i = 5i tN 5j tN - 6^5 

Writing equation 2.33 in indicial nomenclature leads to qtj + q Jt i = 5i t N^j,N ~ l which is the 
desired result. □ 

2.3.1 Collocation 

A Legendre collocation operator may be constructed by approximating the integrals in equations 
2.30, 2.31 and 2.4 by the LGL quadrature formula. Let r] = (rji, 772 , ■ • • t/jv-i, t]n) be the nodes 
of the LGL quadrature formula (i.e., the zeroes of the polynomial P n _ 1 (x)( 1 — x 2 ) [31]), and let 
uji, 1 < l < N be the quadrature weights. 

Next, define new mass and stiffness matrices V and Q c by the expressions 

V = Ez L ( 7 7i; x )[L(r7z;x)] T 77i dx ; Q c = E; L (Vh X )[ L '(^; x)f ry dx . (2.34) 
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The matrix V is SPD for any x. The proof is structurally equivalent to the Galerkin proof given 
in equation 2.31 and is omitted. 

Note that in general, V ^ V. The LGL formula is exact for polynomials of degree 2p — 1 but 
L(x)[L(x)] T dx is of degree 2 p. Thus, the integration differs for the highest order term (i.e. , 
2 p th ). Indeed, the two matrix norms differ by a rank one perturbation, i.e V = V + 'y p 'D p eo['D p eo] 1 
where eo = [1, 0, • • • , 0] T , V p is the highest derivative supported by the polynomial, and depends 
on polynomial order. 

The matrices Q and Q c are equivalent. This follows from the fact of the two matrices are 
defined by the polynomials Jj 1 1 L(x)[L/(x)] dx that have a combined rank of 2p — 1. Therefore, 
integration is exact when using the LGL integration formula. 

The uniqueness of the differentiation matrix T> yields the expression 

V = V~ X Q. = V~ X Q. 

This statement does not contradict the fact that V / V. Indeed, expanding the V using the 
Sherman- Morrison formula demonstrates that the difference in the inverses lies in the null space of 
the singular V matrix. 

2.3.2 Diagonal Norm SBP Operators 

Theorem 2.5. The matrix V is diagonal for collocations points located at the LGL quadrature 
points, i.e., x = rj. Furthermore the diagonal coefficients ofV are the integration weights uii, 1 < 
l < N used in the quadrature. 

Proof. Recall that the Lagrange polynomials evaluated at the knot points satisfy the property 
Li{xj) = Sij. Thus, the result follows immediately from the definition of the norm V = 
Ez L (^;x)[L(?7/;x)] T 77;. □ 

2.4 SAT Penalty Boundary and Interface Conditions 

The imposition of boundary and interface conditions is of critical importance in all numerical meth- 
ods for partial differential equations on finite domains. The manner in which these conditions are 
imposed greatly affects the stability and accuracy of solutions. Accurate, stable, and conserva- 
tive interface coupling techniques are essential in multi-domain settings and become critical as the 
domains (i.e., elements) are refined in size. 2 

A straightforward method that permits formal analysis and maintains design-order accuracy is 
the SAT penalty method [26-30] . The approximate spatial integration of the semi-discretization in 
2 . 20 , 

^l T Pq = /o - In + f { N - fV + lT (Sb + Si ) » ( 2 -35) 

illustrates the purpose of the penalty, that may be thought of as replacing some of the computed data 
in the approximation with known data from the boundary condition to specify a mathematically 
well-posed problem. 

^Indeed, the ratio of interior to interface terms remains constant independent of resolution in spectral element 
discretizations; the ratio goes to zero in a multi-block finite-difference discretization. 
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3 Entropy Stable Spectral Collocation: Single Domain 

3.1 Continuous Analysis 

3.1.1 Smooth Solutions 

Consider a nonlinear system of equations (e.g., the Navier-Stokes equations given in 2.1) and assume 
that the solution is smooth for all time. The objective is to bound the solution as sharply as possible. 
A quadratic or otherwise convex extension of the original equations is sought (the existence is in 
general not guaranteed), that when integrated over the domain depends only on boundary data 
and dissipative terms. Fortunately, the Navier-Stokes equations have a convex extension, referred 
to as the entropy function that provides a mechanism for proving stability of the nonlinear system. 

Definition 3.1. A scalar function S = S(q ) is an entropy function of equation 2.1 if it satisfies 
the following conditions: 

• The function S(q) is convex and when differentiated, simultaneously contracts all spatial fluxes 
as follows 

S q f Xi = S q f q q Xi = F\q Xi =4 ; i = (3.1) 

for each spatial coordinate, d. The components of the contracting vector, S q , are the entropy 
variables denoted as w T = S q . F l (q) are the entropy fluxes in the i-direction. 

• The entropy variables, w, symmetrize equation 2.1 if w assumes the role of a new dependent 
variable (i.e., q = q(w)). Expressing equation 2.1 in terms of w is 

Qt + (f l ) Xi ~ (f (v)l )xi = qwWt + {P w )w Xi - (cijW Xj ) x , = o ; i = l,---,d (3.2) 

with the symmetry conditions: q w = [q w ] T , ff = ffj ,dij = c J], 

Because the entropy is convex, the Hessian S qq = w q is symmetric positive definite, 

C T S qq ( >0, VC ^ 0, (3.3) 

and yields a one-to-one mapping from conservation variables, q, to entropy variables, w T = S q . 
Likewise, w q is SPD because q w = w q ~ l and SPD matrices are invertible. The entropy and corre- 
sponding entropy flux are often denoted an entropy-entropy flux pair, (S, F ) . Likewise, the poten- 
tial and the corresponding potential flux (defined next) are denoted a potential-potential flux pair, 

[U]- 

The symmetry of the matrices q w and ff, indicates that the conservation variables, q, and 
fluxes, f l , are Jacobians of scalar functions with respect to the entropy variables, 

q T = tfiw, iff = (3-4) 

where the nonlinear function, ip, is called the potential and ifl 1 are called the potential fluxes [11]. 
Just as the entropy function is convex with respect to the conservative variables (S qq is positive 
definite), the potential function is convex with respect to the entropy variables. 

The two elements of Definition 3.1 are closely related, as is shown by Godunov [32] and Mock 
[33]. Godunov proves that: 
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Theorem 3.1. If equation 2.1 can be symmetrized by introducing new variables w, and q is a 
convex function of<p, then an entropy function S(q) is given by 


ip = w T q — S, 

and the entropy fluxes F l {q ) satisfy 

if 1 = w T f - F\ 


(3.5) 


(3.6) 


Mock proves the converse to be true: 

Theorem 3.2. If S(q) is an entropy function of equation 2.1; then w 1 = S q symmetrizes the 
equation. 

See reference [34] for a detailed summary of both proofs. 

Entropy analysis is now applied to the Navier-Stokes equations to determine the limits of non- 
linear stability. 

Contracting equation 2.1 with the entropy variables results in the differential form of the entropy 
equation, 


S q qt + S q f(q) Xi = S t + F Xi = S q f £> 






(3.7) 


Integrating equation 3.7 over the domain yields a global conservation statement for the entropy 
in the domain 

d 

w" 


it J Si « - 

n 


■ T f( v ) _ p 


J an 


~ / w x CijW Xi dxi. 


(3.8) 


n 


It is shown elsewhere [9, 10] that dij in the last term in the integral are positive semi-definite. Note 
that the entropy can only increase in the domain based on data that convects and diffuses through 
the boundaries. The sign of the entropy change from viscous dissipation is always negative. 

Thus, the entropy equation derived in 3.8 is the convex extension of the original Navier-Stokes 
equations, and the entropy function serves as an estimator of stability of the system. 


3.1.2 Discontinuous Solutions 


The Euler terms in equation 2.1 (the convective terms to the left of the equal sign) admit discon- 
tinuous solutions in finite time even for smooth initial and boundary data. Thus, weak solutions 
to the integral form of equation 2.1 are appropriate for these situations. Although equation 3.8 is 
an integral statement of entropy conservation, it is not strictly valid in the presence of discontinu- 
ities, because it does not accurately account for the dissipation of entropy at the discontinuity (i.e., 
shocks). 3 Although the precise amount of entropy dissipated at a shock is not known apriori, what 
is known is the sign of the jump in entropy. Thus, a general (yet somewhat ambiguous) statement 
of the conservation of entropy in the domain is 


it J SiXi - 

n 


w T f (v) _ F 


J an 


J 


wlfijWxi 


d Xi 


(3.9) 


3 Note that the mathematical entropy has the opposite sign from thermodynamic entropy in gas dynamics. 
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Weak solutions in general may not be unique [24,35]. In these cases, Equation 3.9 is used to identify 
spurious solutions that violate the entropy condition from those that are physically admissible. 
Entropy analysis is valid for nonlinear equations and even those that admit discontinuous solutions; 
it is more generally applicable than linear energy analysis and gives a stronger stability estimate. 

3.2 Semi-Discrete Entropy Analysis 

The semi-discrete entropy estimate is achieved by mimicking term by term the continuous estimate 
given in equation 3.8. The nonlinear analysis begins by contracting the entropy variables, w T , 
with the semi-discrete equation 2.20. (For clarity of presentation, but without loss of generality, 
the derivation is simplified to one spatial dimension. Tensor product algebra allows the results to 
extended directly to three-dimensions.) The resulting global equation that governs the semi-discrete 
decay of entropy is given by 


w T Pq t + w T Af = w 1 Af 1 ^ + w 1 g;,, (3.10) 

where 

w = (w(q 1 ) T ,w(q 2 ) T , ■ • • ,w(q N ) T ) , 

the vector of entropy variables. Each semi-discrete term is now analyzed to demonstrate that it 
mimics the corresponding term in continuous entropy estimate, provided that a diagonal norm SBP 
operator is used. The form of the penalty terms is presented in a later section 4. 

3.2.1 Time Derivative 

The time derivative is by definition in mimetic form for diagonal norm SBP operators. Arbitrary 
diagonal matrices commute, so the pointwise definition of entropy 

wf{qi)t = (Si) t , Vi. 

yields the expression 

w T Vq t = l T Pw r q ( = l T VS u q t = l T VS t . 

3.2.2 Inviscid Flux Conditions 

The inviscid portion of equation 3.10 is entropy conservative if it satisfies 

w T A{ = F(q N )-F{ qi ) = l T AF. (3.11) 

Recall that w and f , F are defined at the solution points and flux points, respectively. One plau- 
sible solution to equation 3.11 is a pointwise relation between solution and flux-point data, which 
telescopes across the domain and produces the entropy fluxes at the boundaries. Tadmor [11] devel- 
oped such a solution based on second-order centered operators. Herein, this solution is generalized 
for Legendre spectral collocation operators of any order. 

Theorem 3.3. The local conditions 

(rn + 1 - m) T fi = 4>i+i -4>i, i = 1, 2, . . . ,N - 1 ; V’l = V’i> V ’ n = V ’ N (3.12) 
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when summed, telescope across the domain and satisfy the entropy conservative condition given in 

equation 3.11. The potentials ifi+i and ifi need not be the pointwise ipi + 1 and ifi, respectively. A 

—(S) 

flux that satisfies this condition given in equation 3.12 is denoted f . 

Proof. Substituting the definition for generalized summation-by-parts in Section 2.1.4, A = B — A, 
into the global entropy conservation condition in equation 3.11 yields 

w T Rf - w T Af - 1 t BF + 1 t AF = w r Rf - 1 t BF - w T Af = 0. (3.13) 

The boundary terms in 3.13 may be reorganized as 

w r 6f - 1 T B~F = ( wJff N - F n ) - (vjjfi - Fi) = if N - ifi = if 7 ' Bl, (3-14) 

where if\ and i/’jv represent the potential flux defined in equation 3.6. Defining [f] as a diagonal 
(N + 1) X ( N + 1) matrix containing the elements of f, and substituting equations 3.13 and 3.14 
into 3.11 yields 

(if t b - w T A[f]) 1 = 0. 

~ T ~ T ~ - 

Substituting the equality if B1 = if A1 into the left hand side of the equation yields 

(if T A — w T A[f]^ 1 = 0. (3.15) 

This is satisfied by the vector sufficient condition, 

if T A = w r A[f], ifi = ifi, ifN = i’N ■ (3.16) 

A pointwise examination of the vector condition yields the desired result. □ 

Note that Tadmor arrives at the condition if = if because of an assumption on the form of 
F [11]. In the generalized condition derived in equation 3.16 it is unnecessary to define if in the 
domain interior. This generality is important because the entropy conservative fluxes needed for 
high-order methods do not satisfy Tadmor’s original two point form [12]. 

The local conditions defined in equation 3.12 are sufficient for three-point, second-order centered 
discretizations. The entropy conservative flux is constructed at each flux point, /,;, from solution 
point data obtained from the two adjacent solution points, i and i + 1. The resulting operator 
telescopes across the domain when contracted with the entropy variables. 

Likewise, equation 3.12 is generalizable to coarser grids using the notion of Richardson extrap- 
olation; e.g., two point, second order, skew-symmetric operators spanning 5 pt, 7 pt, 9 pt, • • • . Thus, 
if a high-order derivative operator may be constructed from linear combinations of these elemen- 
tal second-order operators, then high-order accuracy is achieved that ensures entropy consistency. 
(Schemes of this class are devised for periodic domains by LeFloch and Rhode [12], and extended 
to high-order periodic domains by LeFloch and Rhode [12].) 

The matrix footprint of any spectral collocation operator is dense and involves (p+ l) 2 indi- 
vidual terms. Neither of the two previous approaches extend in general to these dense operators. 
A different strategy for constructing high-order entropy conservative fluxes is presented in refer- 
ence [9, 10], and utilizes linear combinations of two-point entropy conservative fluxes, combined 
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using the coefficients in the SBP matrix Q. This new approach follows immediately from the gen- 
eralized telescoping structural properties of diagonal norm SBP operators given in section 2.1.4. 
Because is requires only the existence of a two-point entropy conservative flux formula and the co- 
efficients of the Q, it is valid for any SBP operator which satisfies the constraints given in equation 
2.3. Thus, it is valid for Legendre spectral collocation operators. 

The proofs of this alternative approach for building entropy conservative operators of any order 
are quite involved. For brevity only two of the theorems are included herein. Interested readers 
should consult references [9, 10] for details. 

The first theorem establishes the accuracy of the new fluxes: that a high-order flux constructed 
from a linear combination of two-point entropy conservative fluxes retains the design order of the 
original discrete operator for any diagonal norm SBP matrix Q. 

Theorem 3.4. A two-point entropy conservative flux can be extended to high order with formal 
boundary closures by using the form 

N i 

f(.S) = Y Y 2 q(e,k)fs (qe, qk), 1 < i < N - 1, (3.17) 

k=i-\- 1 £=1 

when the two-point non- dissipative function from Tadmor [11] is used 

i 

Is (qk, qe) = J g ( w(q k ) + £ (w(q e ) - w(q k ))) d£, g(w(u )) = f(u). (3.18) 

o 

The coefficient, qi k ,e)> corresponds to the (k,£) row and column in Q, respectively. 

Proof. To show the accuracy of approximation, the flux difference is expressed as 

N i N 2—1 

ip-m= £ e 2 q(e,k)fs (ug, u k ) ££ 2 q(e,k)fs (ue, u k ) , 2 < i < N — 1. 

k=i+l £=1 k=i £=1 

that may be manipulated into the form (see references [9,10]) 

N 

fi S) ~ fill = 2 q(iJ)fs(ui,Uj), 1 <i<N. (3.19) 

3 = 1 

This form facilitates an analysis by Taylor series at every solution point by using the expression 
for the two-point fluxes given in equation 3.18. The remainder of the proof is presented elsewhere 
[9,10], □ 

The second theorem establishes that the linear combination does indeed preserve the property 
of entropy stability for any arbitrary diagonal norm SBP matrix Q. 

Theorem 3.5. A two-point high-order entropy conservative flux satisfying equation 3.12 with for- 
mal boundary closures can be constructed using equation 3.17, 

N i 

f!' S) = Y 2( ip,k)fs (qe, qk), i < i < N - l, 

k=i-\- 1 £=1 
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where fs{Qi,Qk ) is any two-point non-dissipative function that satisfies the entropy conservation 
condition 

{w t - w k ) T Js ( qe , Qk) = ipe - fi>k- (3.20) 

The high- order entropy conservative flux satisfies an additional local entropy conservation property, 

w r ?-'Af (S) = V-'AF = F x ( q) + T d , (3.21) 

or equivalently, 

(f? ] - ftl) = (Fi - ^-i) , 1<*<V, 

where 

N i 

Fi = ^(4*0 + ™fc) T /s (®> 9fe) - (# + V 7 *:)] , 1 < i < iV - 1 

/c=z+l ^=1 

Proof. For brevity, the proof is not included herein, but is reported elsewhere [9,10]. □ 

.Remark. The existence of a local second-order entropy flux satisfying the two point shuffle 
relation given in equation 3.20 is a very strong constraint, and has until recently been a computation 
bottleneck [13], 

Remark. The entropy consistency proof is satisfied for all two-point fluxes that satisfy 3.20. 
The accuracy proof is proven only for fluxes in the integral form 3.18. Currently, the proof does 
not extend to any flux satisfying 3.20, so such fluxes should be validated for accuracy independent 
of Theorem 3.4. 

In summary, Theorems 3.4 and 3.5, guarantee that the extension of the two-point flux given in 
equation 3.17, is a high-order accurate entropy conservative discretization of the conservation law. 

3.3 Entropy Stability of the Euler Equations 

Tadmor [11] identifies three “tools of the trade” in the analysis of entropy stability: comparison ar- 
guments, a homotropy approach, and kinetic formulations. In a comparison approach, the entropy 
dissipation generated by the primary scheme is compared with the baseline entropy of a scheme 
known to be at least entropy conservative. If the dissipation is less than the entropy conservative 
datum, then more dissipation is necessary. Conditions that guarantee entropy stability are devel- 
oped in references [9, 10], and extend directly to this work. Stabilization of the inviscid terms is 
the topic of ongoing work. 

3.4 Entropy Stable Viscous Terms 

Using the formalism introduced in Section 2.1.2, viscous terms are derived that discretely mimics 
the continuous entropy properties. As with the continuous estimate, the proof requires the viscous 
fluxes to be written as functions of the discrete gradients of the entropy variables, 

(caw x ) x = 'P - 1 Af^ = V 2 {cii ) w + T p , ^ 

V 2 (c ii )w = V- 1 {-V T V[c ii }V + B[c ii ]V) w . 


(3.22) 

(3.23) 
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The accuracy requirements are automatically satisfied. The coefficient matrix [c&] is positive semi- 
definite because it is constructed from block-diagonal combinations of positive semi-definite matri- 
ces. 

The contribution of the viscous terms to the semi-discrete entropy decay rate is 

w T Af^ = w,8[cii]Z>w — (Vw) T 'P[cn\(Vw). (3.25) 

The last term is negative semi-definite. As with the continuous estimate given in 3.8, only the 
boundary term can produce a growth of the entropy, and thus the approximation of the viscous 
terms is entropy stable. (Wellposed boundary conditions bound these terms.) 

The high-order entropy stable viscous terms described in equations 3.24 and 3.25 first appear in 
a companion work [9, 10], that focuses primarily on the entropy stability of WENO finite-difference 
operators. Unlike multi-domain finite-difference formulations, spectral element operators have the 
potential to superconverging from p to p + 1, with the appropriate choice of viscous interface 
coupling terms. Thus, a slightly modified discretization of the viscous terms is adopted herein. 


4 Entropy Stable Spectral Collocation: Multi-Domains and Dis- 
continuous Interfaces 

Design order consistent, conservative interface treatments that are provably entropy stable are de- 
veloped next. Two popular techniques for coupling interfaces are the discontinuous and continuous 
approaches. The focus herein is on the discontinuous approach. 

Numerous discontinuous coupling approaches exist for spectral element operators. See refer- 
ence [36] for a review of the (weak form) finite element literature or [18] for a general treatment 
of strong form collocation approaches. The advantage of discontinuous approaches is that each 
domain may be treated individually. The inviscid terms are typically coupled across the interface 
through a Riemann solver, a technique that ensures conservation and which provides a necessary 
mechanism for dissipating unresolved modes in the simulation. Likewise, numerous approaches 
exist for coupling the viscous terms across interfaces, but two of the more popular are the Local 
Discontinuous Galerkin (LDG) FEM approach proposed by Cockburn and Shu [37], and the second 
technique proposed by Bassi and Rebay [38]. The present work develops a new entropy stable LDG 
approach based on the strong (differential) form of the governing equations. 


4.1 Navier-Stokes in One Spatial Dimension 

For clarity of presentation we analyze the one-dimensional case. Extension to three spatial di- 
mensions is straightforward although algebraically involved. The inviscid and viscous penalties 
are developed separately to guarantee an entropy stable inviscid penalty in the limit of vanishing 
viscosity. 

A spectral collocation approximation for the Navier-Stokes equations, including interface cou- 
pling terms, motivated by the LDG FEM approach is 
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(4.1) 

with the subscripts l,r denoting the “left, right” elements (subdomains). The subscript i denotes 
an interface quantity with the superscripts “(-),(+)” denoting the collocated values on the left and 
right side of the interface, respectively. The flux i ssr (q^\ g( + )) is an entropy stable reconstructed 
flux. 


Theorem 4.1. The approximation (4-1) of the one- dimensional Navier-Stokes equations is entropy 
stable if the reconstructed flux f ssr (q ( -~\ q (+)) and the viscous matrix parameters Loo,Loi,L\o, 
Roo,Rqi,Rio satisfy the following conditions. The entropy stable reconstruction f ssr (q(~\ gW) must 
be more dissipative than the entropy conservative reconstruction. A sufficient condition is the use 
of dissipation of the form 

- w (_) ) f ssr (g(~) ,£(+)) = ^(+) - -i/iH -T diss . (4.2) 

where T d i ss is an interface term which is uniformly dissipative (e.g., Lax- Friedrichs dissipation), or 
zero. 

The viscous parameters must satisfy the conditions. 


Loo = Roo < 0 ; iioi = +Toi + c u ; Rio = L io + c n ; L\o = — Toi — c u ; 

(4.3) 

where symmetric dissipation matrix cu is defined as 

« = Ki"’ + i +) ) ■ ( 4 - 4 ) 

Finally, suitable boundary conditions and initial data must be provided. 


Proof. Entropy stability of (4.1) follows if the interface treatment at x = Xi is more dissipative 
than the entropy conservative interface treatment. The farfield boundary penalties are assumed to 
be stable, allowing separate analysis of the interface terms. 

The entropy method is used to prove the stability of equation (4.1). Multiplying the two discrete 
equations in the left subdomain by wf and ( cm Ei) t , respectively, and the two discrete equations 
in the right subdomain by wj and ( ca r E r ) T , respectively, and then summing the four equations 
and collecting terms results in the expression 


where 




,2 

\Vr 


= T 


(4.5) 


T i = [2 wiCiu Ei - F(q)f_ 1 + [2 w r c iir E r - F(q)}} + 

+ ~ f ssr (<? (_ W +) )) ~w\ +) {F(ql +) ) - f ssr (g (_) , g (+) )) 
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(4.6) 
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The viscous dissipation terms || y/cm H ; \\^ and \\\fcn r S r ||^ are uniformly dissipative. Thus, en- 
tropy stability of equation (4.5) follows immediately if the term Tj is dissipative. Note that Tj is 
composed of both inviscid and viscous terms; i.e. , T j = 'Y I . nmsid + 'fViscous . The inviscid and 
viscous terms are bounded individually to guarantee that the inviscid terms are stable in the limit 
of Re — > oo. 

Inviscid Stability 

The inviscid interface terms in Tj are 

X Invisid _ _ F(q ( '~' ) ) 

+ - f Mr (?C-W+))) - ^ +) (F(^ +) ) - P-(gH, ?(+>)) (4 ' 7) 

Substituting the definitions for the entropy fluxes 

F{q- +) ) = w\ +) f(q\ +) ) + ^ (+) ; R(gJ _) ) = w\~ ] /(gj _) ) + ^ , 

into equation 4.7 yields the equation 

T{ nvisid = (w^ -w^Y f ssr {q(-\q^) - (V> (+) - V’ (_) ) (4.8) 


that is a dissipative term provided that condition 4.2 is met. 

Viscous Stability 

Stability of the viscous interface terms is proven by collecting all remaining terms in 4.6 and 
expressing them as a matrix contraction T Y lscous = jF Mi % , where 
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Mi 


(+2Loo) 

— {L 00 + i?oo) 

( cm + Toi + Too) 

— (Lqi + Rio) ( 


— (Lqo + T?oo) 
(+2i?oo) 

— (Rqi + Lio) 
Ciir + Roi + T?io) 


(cm + Toi + L io) 
— (T?oi + Lio) 

(2cqCjj/) 

0 


— (Lq\ + R\o) 

( Ciir + T?oi + Rio) 
0 

(2 Oi r Ca r ) 


(l- 9 ) 

Note that the matrix Mj is symmetric by construction. Stability of the interface follows immediately 
if Mj is negative semi-definite, a matrix property that is easily established by testing for the 

existence of a Cholesky factorization. The admissible values of the six viscous parameters Loo -Rio 

are determined by requiring strictly negative diagonal pivots during the Cholesky factorization of 
Mi which yields the desired result. □ 


Conservation is a desirable property of any numerical method based on the strong form of the 
equations, and guarantees that convergent solutions satisfy the weak form of the governing equations 
away from discontinuities. The conservation property is established by multiplying equation 4.1 by 
an arbitrary smooth test function, and then transferring the action of the derivative operator (via 
matrix manipulations) over to the test function. 

All inviscid terms at the interface vanish immediately by the specific choice of the penalty, 
provided that the test function is continuous across the interface. The viscous condition 


Lq\ = Rqi — Cjj 
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is necessary for the approximation (4.1) to retain conservation (in the V norm) across the interface. 
Note that this condition is also identified as a necessary constraint during the Cholesky factorization 
of the matrix M % . 

4.1.1 Inviscid interface dissipation 

The entropy stable flux f ssr (gl~). q (+)) used in equation 3.12 is composed of the sum of an entropy 
conservative flux i sr (q^~\q^) and a dissipation term . Lax-Friedrichs dissipation is perhaps 
the simplest means of stabilizing the interface. The resulting flux may be expressed as 

f ssr (g( - ), g( + )) = f sr (q(~\ </(+)) + — w^) , 

A = [i((u(-)) 4 + (c(+)) 4 + («(-)) 4 + (c«) 4 )]* . (4 ' 10) 

A flux that includes dissipation in this form is denoted an entropy stable Lax-Friedrichs flux. Note 
that this form differs from the conventional Lax-Friedrichs flux by replacing the linear average 
of the interface fluxes \{f{q^) + f(q ^)), with the nonlinear average entropy conservative flux 
f sr (<?(-) ! (/("I - )). In practice, the dissipation term may be scaled by any positive number. The 
resulting term produces in the entropy estimate as the interface terms are contracted against the 
entropy variables is of the form 

— A(u/ - ) — w^) 2 (4-11) 

Lax-Friedrichs dissipation overly damps the convective waves at the interface. A more refined 
approach to dissipate scales each characteristic component based on the magnitude of its eigenvalue. 
A flux that includes dissipation of this form is denoted an entropy stable characteristic flux, and is 
implemented as 


f aar ( 9 H , 9 (+) ) = f sr (q ( ~\q {+) ) + l/2y\\\y T {w^ - «,«) , 

(f) q = yw T ; = yy T [ ' J 

Note that the relation q w = yy T is achieved by an appropriate scaling of the rotation eigenvectors. 
See the work of Merriam [39] for more details. 


5 Entropy analysis for the Navier-Stokes equations 

5.1 Euler and Navier-Stokes Equations 

The end goal of this work is to simulate the Navier-Stokes equations in a stable and efficient manner. 
Toward this end, the application of entropy stable inviscid and viscous terms are detailed in this 
section. The conservative variables for the Navier-Stokes equations are 

q = (P; pv \ , pv 2 ,pv :h f)E) T , (5.1) 

where p denotes density, v = {v\,V 2 ,V 3 ) T is the velocity vector, and E is the specific total energy. 
The convective fluxes are 

f = ( pvi , pviv i + Sap , pviv 2 + Si2P, pviv 3 + S i3 p, pviH) T , (5.2) 
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where p represents pressure, H = E + p/p is the specific total enthalpy, and 5ij is the Kronecker 
delta. The viscous flux terms are 


/b’b = (0, Tii,T i2 , Tj3, TjiVj - qi) T , 

(5.3) 

where the shear stress is 


T ij = T f ( v i)xj T {vj)xi — j y 

(5.4) 

and the heat flux is 

Qi — . 

(5.5) 

T denotes the static temperature, with p = p(T) and k = k(T) the dynamic viscosity and thermal 
conductivity, respectively. 


The viscous terms may also be expressed as 


f (v)l = Cijq Xj , (5.6) 

which is a convenient form for entropy analysis. The constitutive relations for a perfect gas are 

h = H — -VjVj = c p T, (5-7) 

where c p is the constant specific heat, and 

v=pRT ' R= ~m’ (5 8) 

where R u is the universal gas constant and MW is the molecular weight of the gas. The speed of 
sound for a perfect gas is 

c=^RT, 7 = T“lv (5-9) 

Cp it 

In the entropy analysis that follows, the definition of the thermodynamic entropy is the explicit 
form, 

s = ^T log (0- filog («)' <510) 

where To and po are the reference temperature and density, respectively. 

5.1.1 Entropy Analysis 

In the Navier-Stokes equations, entropy stability is not the same as full non-linear stability. Nev- 
ertheless, entropy stability gives a stronger stability estimate than linear energy stability, and in 
many ways is easier to apply. In this section, the continuous entropy stability analysis is conducted 
first to illustrate the entropy characteristics of the governing equations. Discrete spatial operators 
are derived next, via semi-discrete entropy analysis, that mimic these continuous properties. 

The entropy-entropy flux pair for the Navier-Stokes equations is 

S = —ps, F l = —pviS, (5-11) 


22 



and the potential-potential flux pair is 


f = pR, if) 1 = pviR. 


(5.12) 


Note again that the mathematical entropy has the opposite sign of the thermodynamic entropy. 
To avoid confusion, herein entropy refers to the mathematical entropy unless otherwise noted. The 
entropy variables using the pair in 5.11 are 


w = 


St 



2 T 


V\ V 2 

) rj~i '■ 


V3 

T ’ 



(5.13) 


and may be shown to have a one-to-one mapping with the conservative variables provided p,T > 0. 
Expressly: 

C T SggC T > 0, VC + 0, p,T> 0. 

This restriction is what makes the entropy proof fail to be a true measure of nonlinear stability. 
Another mechanism must be employed to bound p and T away from zero to guarantee stability. 

The entropy equation is found by premultiplying the Navier-Stokes equations with the transpose 
of the entropy variables, 


SqQt + S q (f l ) Xi = S t + F*. = w T ( Cijq w w Xj ) Xi = ( w T Cijq w w Xj ) Xi - w x .Cijq w w Xj , 
where the viscous terms satisfy 

C-ij^w = Cij = Cjj, Cxi^ijCxj ^ 0, VC, 

provided 

T > 0, p(T) > 0, k(T) > 0. 

The total entropy decay rate is found by integrating 5.14 over space, 


(5.14) 


(5.15) 


j t JsdV = J {w T c ij w Xj - F i ) dS i - j 


an 


- I w x .CijW Xj dV. 

n 


(5.16) 


5.1.2 Discretization Notes 

To facilitate the extension of the entropy stable methods to the three-dimensional equations, we 
define the three-dimensional nomenclature and examine the general form of the semi-discretization. 
The semi-discrete form of 2.1 is 

q* + E ( F - = E (st + g/) • ( 5 - 17 ) 

i = 1 2—1 

The solution vector is ordered as 

q = ('U(X( 1 )( 1 )( 1 )) T , u{x( 1)(1)(2)) T , ■ • ■ , u{x {Nl)(N2)i N 3 )) T ) T ■ 

The roman superscript indices on the flux in 5.17 indicate the direction of the flux, and parenthetic 
superscripts indicate the type of flux, i.e. V for viscous and S for entropy conservative. 
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5.1.3 Entropy Stable Spatial Discretization 


The inviscid terms in the discretization of the Navier-Stokes equations are calculated according to 
equations 3.17 and 3.18, by using the two-point entropy conservative flux of Ismail and Roe [13], 

Qi+ 1 ) = {pVj,pVjVi + Sjip, pVjV 2 + Sj2P, pVjV 3 + Sj 3 p, pvjH)^j , 


Vj I Vj+l 

V% ^ 


v = 


Pi I Pi+i 
V77 y/Ti + l 


1 _|_ 1 ’ P — _J 1 ’ 


VT 


log 


h = R- 


VTi.pi 


y/Ti+iPi+l 


i + 

VTi ^ 


VT 

y/TiPi + y/Ti+ipi+i 


v ( A + vfc) ^ Tlpl ~ (5-18) 


7 + 1 


log ( \/W 


7 - 1 


log 


Tj Pi 


H = h + -vivi, p = 


^4+1 Pi+1 J \^y/Ti yjTi- |_i J ) 

+ 7Cr) ('fTiPi - V T i+\Pi+l) 

2 (\og(y/TiPi) ~ \og{y/T i+ ip i+1 )) 


This somewhat complicated explicit form is the first entropy conservative flux for the convective 
terms with low enough computational cost to be implemented in a practical simulation code. Pre- 
viously, Tadmor [2] derived an entropy conservative flux form that required integration through 
phase space, but this was deemed too expensive to be practical. 


5.1.4 Energy Stable Boundary Conditions 

The problems studied herein are limited to open boundary conditions, which are implemented as 
suggested by Svard et al. [28] 


6 Accuracy Validation 

The new entropy stable methodology is validated using linear advection, nonlinear Burgers’ equa- 
tion, the Euler equations, and the Navier-Stokes equations. The methodology is also extended to 
three dimensions for the compressible Euler and Navier-Stokes equations, the target applications 
for this work. 

6.1 Test Equations 

The accuracy and robustness of the algorithms developed herein are tested using two smooth and 
two discontinuous problems. The smooth problems are the propagation of an isentropic vortex, and 
the propagation of the viscous shock. Both problems demonstrate the design-order convergence of 
the new entropy conservative formulation. 
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6.1.1 Linear Advection 


Linear advection tests the accuracy of the first derivative terms. Design order convergence of p+ 1 
is sought in all cases. The linear advection equation is given by 


du | dc u 
dt ' dx 


= 0 ; 


-1 < x < 1 ; t >= 0 , 
= 9-1 (t) , 


( 6 . 1 ) 


with an exact solution given by 

u(x,t ) = exp(— — ct ) 2 ) ; c = 1 . ( 6 . 2 ) 

Initial and boundary data is provided that is consistent with the exact solution. 


6.1.2 Burgers’ Equation 

The nonlinear Burgers’ equation is a one-dimensional model for the inviscid-viscous interaction 
found in the full Navier-Stokes equation. Burgers’ equation is given by 


du _i_ 1 du 2 
dt ' 2 dx 

au(—l,t) — €U x (—l,t) 


€0 ; —l<x<l;t >= 0 , 

9-i(t) , (3u(-l,t) + eu x (-l,t) = g-i(t) . 


An exact solution is given by 


v(r t) - acx P(( b - a )^ - ct - d )/(^)+ b . „ - _I u - i r - i( n 4 . h) d - i 

“l-DD - lexp((6-o)(x-ct-d)/(2e)+l ’ 2 , V — ±, C — 2 {U U) , U — 2 

Initial and boundary data are provided that is consistent with the exact solution. 


(6.3) 


(6.4) 


6.1.3 Isentropic Vortex 

The isentropic vortex is an exact solution to the Euler equations and is an excellent test of the 
accuracy and functionality of the inviscid components of a Navier-Stokes solver. It is fully described 
by 

f(x,y,z,t ) = 1 - {x - x 0 - Uoo cos(a) t) 2 + (y - y 0 - sin (a) t) 2 , 

T(x,y,z,t) = 1 - elM^^e^p(f(x,y,z,t)) , p(x,y,z,t ) = TV i, 

u(x,y,z,t ) = Uoo cos(a) - e v y ~ y °~ U £ sm{a)t exp \ ; 

v(x,y,z,t ) = Uoo sin(a) - e^ x ~' ro ~^ cos(a)t exp , 

w(x,y,z,t) = 0. 

In this study the values U 0 0 = MooCoo , e v = 5.0, = 0.5, and 7 = 1.4 are used. 

The Cartesian grid test case is described by 


(6.5) 


x € (—15, 15), y € (—15, 15), (x 0 ,y 0 ) = (0,0), a = 0.0, t > 0. 
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6.1.4 The Viscous Shock 


The Navier-Stokes equations support an exact solution for the viscous shock profile, under the 
assumption that the Prandtl number is Pr = Mass and total enthalpy are constant across a 
shock. Furthermore, if Pr = ^ then the momentum and energy equations are redundant. The 
single momentum equation across the shock is given by 


avv x — (v — l)(u — Vf) = 0 


— OO < X < oo 

v = jl. ■ Vf = m ■ a = 7=1 v 
ul, ’ U J ill ’ 27 Pr m 


t > 0; 


An exact solution is obtained by solving the momentum equation for the velocity profile. 

h a (■ L °g\( v - l)(v - v f )\ + ^Log\^\) 


X = 2 C 


( 6 . 6 ) 


(6.7) 


A moving shock is recovered by applying a uniform translation to the solution. A full derivation 
of this solution appears in the thesis of Fisher [22], 


6.2 Test Results 

All tests include elements of polynomial degrees 1 < p < 4. A uniform grid refinement study is 
performed using a grid-doubling procedure. When possible, a randomly distributed grid is used. 
A properly nested set of uniformly refined random grids is generated as follows. First, a random 
grid is generated at the coarsest resolution. This grid is then scaled and replicated 2 s , 3 < s < 8 
times on the interval — 1 < x < 1. Thus, the randomness of the coarsest grid is preserved on all 
levels. All simulations use a fourth-order low-storage Runge-Kutta scheme to advance the solution 
in time. A temporal error controller is used to guarantee that the temporal error is subordinate to 
the spatial error, and to automate the determination of the maximum stable timestep for each test 
case. 


6.2.1 Linear Advection 

Table 1 provides data from a uniform grid refinement study of the linear advection equation. Design 
order convergence (i.e., p + 1) is achieved in both the L 2 and L ^ norms. This test verifies that 
the basic spectral collocation first derivative operator super converges by one order, relative to the 
polynomial order of the approximation. 

6.2.2 Nonlinear Burgers’ Equation 

Table 2 contains data from both a uniform and nonuniform grid refinement study of the nonlinear 
Burgers’ equation. The nonuniform grid refinement study is included to identify superconvergence 
resulting from fortuitous cancellation of viscous error terms at element interfaces. (See reference [18] 
for a discussion of this phenomena.) 

Sharp design order convergence of p + 1 is achieved in both the L 2 and Loo norms on the 
uniform grid for polynomials of even order. This sharp convergence may be the result of fortuitous 
cancellation of errors. Design order convergence is achieved asymptotically for polynomials of 
odd order, in both L 2 and L^. The nonuniform refinement study shows asymptotic design order 
convergence for both even and odd polynomial orders. 
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This test verifies that the superconvergence observed in the linear advection study, extends to 
nonlinear inviscid terms, and that the entropy stable LDG implementation is design order accurate 
for the viscous terms. Furthermore, superconvergence may be achieved on irregular grids. 

6.2.3 The Euler Vortex 

The convergence rate for the isentropic Euler vortex is evaluated on a properly nested sequence of 
uniform two-dimensional grids. The vortex profile is initially located in the middle of the domain 
and is simulated until t = 0.25. The reference Mach number is M = 0.5, and the translation 
velocity of the vortex is unity. The errors for the uniform grids are shown in Table 3. 

Theorem 3.4 proves that design order entropy conservative fluxes may be constructed using a 
linear combination of two-point entropy fluxes. A critical assumption used in the proof is that the 
two-point non-dissipative fluxes satisfy Tadmor’s integral relation given in equation 3.18. Herein, 
the non-dissipative Euler fluxes of Ismail and Roe [13] are used. The study provides evidence that 
the non-dissipative Euler fluxes of Ismail and Roe [13] do not degrade the formal accuracy. The 
interfaces are treated by using the entropy stable characteristic fluxes with a weighting parameter 
tuned to produce “upwind fluxes” at the interfaces. Design order convergence is achieved in all 
cases. 

6.2.4 The Viscous Shock 

The convergence rate for the viscous shock is evaluated on a properly nested sequence of uniform 
and nonuniform grids. The shock profile is initially located in the middle of the domain and is 
simulated until t = 1.00. The Reynolds number is Re = 10 and the reference Mach number is 
M = 2.5. The errors for the uniform and the nonuniform grids are shown in Table 4. 

The interfaces are treated by using the entropy stable characteristic fluxes with a weighting 
parameter tuned to produce “upwind fluxes” at the interfaces. Design order convergence is achieved 
in both the L 2 and norms on uniform and nonuniform grids. 

7 Conclusions 

High-order entropy stable spectral collocation methods of arbitrary polynomial order, are derived 
for the nonlinear conservation laws of the Navier-Stokes equations. The discrete operators are for- 
mulated by representing conventional spectral operators in the summation-by-parts framework, for 
which entropy stability of the Navier-Stokes equations has already been established. The individ- 
ual spectral elements are coupled together in an entropy stable and conservative fashion by using 
existing SBP-SAT operators for the convective terms. A new entropy stable local discontinuous 
Galerkin scheme is developed to couple the viscous terms. The new operators are shown to preserve 
design-order accuracy of p + 1 on model problems. 
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Appendix A 


Differentiation Operators 


The discrete operators are provided for polynomial orders one to four in Table Al. Three quan- 
tities completely define the discrete operators; the diagonal norm, V, the nearly skew-symmetric, 
Q and the positions of the collocation points, x. (Recall the the differentiation matrix is defined 
as D = P -1 Q.) The upper triangular portion of Q is only provided. The full Q matrix may be 
reconstructed from the skew-symmetry property Q + Q T = B. 
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Table 1 Error convergence is shown for the linear advection equation. 


P=1 

L z error 

L z rate 

L°° error 

L°° rate 

8 

3.41E-02 


5.09E-02 


16 

7.80E-03 

-1.92 

1.63E-02 

-1.64 

32 

1.99E-03 

-1.97 

4.42E-03 

-1.88 

64 

5.02E-04 

-1.98 

1.14E-03 

-1.95 

128 

1.26E-04 

-1.99 

2.89E-04 

-1.97 

256 

3.16E-05 

-1.99 

7.28E-05 

-1.99 

P=2 

L z error 

L z rate 

L°° error 

L°° rate 

8 

1.85E-03 


2.64E-03 


16 

2.46E-04 

-2.90 

3.46E-04 

-2.93 

32 

3.12E-05 

-2.97 

4.69E-05 

-2.88 

64 

3.91E-06 

-2.99 

5.99E-06 

-2.97 

128 

4.89E-07 

-2.99 

7.56E-07 

-2.98 

256 

6.12E-08 

-2.99 

9.62E-08 

-2.97 

P=3 

L z error 

L z rate 

L°° error 

L°° rate 

8 

5.28E-05 


1.33E-04 


16 

3.09E-06 

-4.09 

9.29E-06 

-3.83 

32 

1.89E-07 

-4.02 

6.04E-07 

-3.94 

64 

1.18E-08 

-4.00 

3.92E-08 

-3.94 

128 

7.36E-10 

-4.00 

2.53E-09 

-3.95 

256 

4.60E-11 

-4.00 

1.59E-10 

-3.99 

P=4 

L z error 

L z rate 

L°° error 

L°° rate 

8 

2.90E-06 


5.62E-06 


16 

9.12E-08 

-4.99 

1.92E-07 

-4.87 

32 

2.83E-09 

-5.00 

6.02E-09 

-4.99 

64 

8.81E-11 

-5.00 

1.86E-10 

-5.01 

128 

2.77E-12 

-4.99 

5.79E-12 

-5.00 

256 

1.17E-13 

-4.56 

5.45E-13 

-3.40 
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Table 2 Error convergence is shown for the nonlinear Burgers’ equation. 




Uniform Grid 



Nonuniform Grid 


P=01 

L z error 

L z rate 

L°° error 

L°° rate 

L z error 

L z rate 

L°° error 

L°° rate 

8 

1.16E-02 


2.40E-02 


1.77E-02 


2.94E-02 


16 

4.43E-03 

-1.38 

9.67E-03 

-1.31 

5.65E-03 

-1.64 

1.25E-02 

-1.24 

32 

1.55E-03 

-1.51 

3.56E-03 

-1.44 

1.87E-03 

-1.59 

4.36E-03 

-1.51 

64 

4.94E-04 

-1.65 

1.19E-03 

-1.57 

5.66E-04 

-1.72 

1.41E-03 

-1.62 

128 

1.46E-04 

-1.75 

3.74E-04 

-1.67 

1.66E-04 

-1.76 

4.56E-04 

-1.62 

256 

4.13E-05 

-1.82 

1.13E-04 

-1.73 

4.69E-05 

-1.82 

1.35E-04 

-1.75 

512 

1.13E-05 

-1.86 

3.29E-05 

-1.77 

1.29E-05 

-1.86 

3.99E-05 

-1.76 

1024 

3.04E-06 

-1.89 

9.36E-06 

-1.81 

3.44E-06 

-1.90 

1.12E-05 

-1.82 

P=02 

8 

1.00E-03 


2.42E-03 


2.04E-03 


3.91E-03 


16 

1.14E-04 

-3.13 

3.08E-04 

-2.97 

2.20E-04 

-3.21 

5.25E-04 

-2.89 

32 

1.42E-05 

-3.01 

3.83E-05 

-3.00 

4.23E-05 

-2.37 

1.16E-04 

-2.17 

64 

1.78E-06 

-2.99 

4.75E-06 

-3.00 

6.54E-06 

-2.69 

1.82E-05 

-2.67 

128 

2.23E-07 

-2.99 

5.94E-07 

-2.99 

1.00E-06 

-2.70 

3.02E-06 

-2.58 

256 

2.79E-08 

-2.99 

7.44E-08 

-2.99 

1.50E-07 

-2.73 

4.89E-07 

-2.62 

512 

3.49E-09 

-2.99 

9.31E-09 

-2.99 

2.19E-08 

-2.78 

7.61E-08 

-2.68 

1024 

4.36E-10 

-2.99 

1.16E-09 

-2.99 

3.07E-09 

-2.83 

1.14E-08 

-2.73 

II 

o 

CO 

8 

7.83E-05 


2.01E-04 


2.23E-04 


5.16E-04 


16 

7.87E-06 

-3.31 

2.41E-05 

-3.05 

2.26E-05 

-3.30 

5.69E-05 

-3.18 

32 

7.55E-07 

-3.38 

2.30E-06 

-3.38 

1.49E-06 

-3.92 

4.47E-06 

-3.66 

64 

6.79E-08 

-3.47 

2.09E-07 

-3.46 

1.26E-07 

-3.56 

4.02E-07 

-3.47 

128 

5.68E-09 

-3.57 

1.76E-08 

-3.57 

1.04E-08 

-3.59 

3.30E-08 

-3.60 

256 

4.43E-10 

-3.67 

1.38E-09 

-3.66 

8.05E-10 

-3.68 

2.62E-09 

-3.65 

512 

3.26E-11 

-3.76 

1.02E-10 

-3.75 

5.91E-11 

-3.76 

1.95E-10 

-3.74 

1024 

2.31E-12 

-3.81 

7.59E-12 

-3.75 

4.14E-12 

-3.83 

1.38E-11 

-3.81 

II 

o 

4^ 

8 

5.51E-06 


1.28E-05 


1.26E-05 


3.53E-05 


16 

1.22E-07 

-5.49 

4.57E-07 

-4.80 

4.72E-07 

-4.74 

1.67E-06 

-4.39 

32 

3.60E-09 

-5.08 

1.32E-08 

-5.11 

3.39E-08 

-3.79 

1.20E-07 

-3.80 

64 

1.11E-10 

-5.01 

4.13E-10 

-5.00 

1.66E-09 

-4.35 

5.83E-09 

-4.36 

128 

3.66E-12 

-4.92 

1.37E-11 

-4.91 

7.69E-11 

-4.42 

2.87E-10 

-4.34 

256 

2.83E-13 

-3.69 

1.67E-12 

-3.03 

3.38E-12 

-4.50 

1.32E-11 

-4.44 

512 

3.08E-14 

-3.20 

1.19E-13 

-3.81 

1.41E-13 

-4.58 

6.00E-13 

-4.45 

1024 

3.65E-14 

0.24 

7.21E-14 

-0.71 

3.67E-14 

-1.94 

7.79E-14 

-2.94 
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Table 3 Error convergence is shown for the isentropic Euler vortex equation. 


Uniform Grid 

P=01 

L z error 

L z rate 

L°° error 

L°° rate 

04x02 

1.42E-02 


5.54E-02 


08x04 

5.87E-03 

-1.27 

2.24E-02 

-1.30 

16x08 

1.88E-03 

-1.64 

7.02E-03 

-1.67 

32x16 

3.98E-04 

-2.24 

1.43E-03 

-2.29 

64x32 

6.38E-05 

-2.64 

2.32E-04 

-2.61 

II 

o 

to 

04x02 

1.83E-03 


1.16E-02 


08x04 

3.74E-04 

-2.28 

2.46E-03 

-2.23 

16x08 

6.40E-05 

-2.54 

4.06E-04 

-2.60 

32x16 

1.20E-05 

-2.41 

7.15E-05 

-2.50 

64x32 

2.69E-06 

-2.15 

1.59E-05 

-2.16 

P=03 

04x02 

1.43E-04 


7.92E-04 


08x04 

1.30E-05 

-3.46 

1.03E-04 

-2.94 

16x08 

9.41E-07 

-3.78 

7.92E-06 

-3.69 

32x16 

6.79E-08 

-3.79 

8.27E-07 

-3.25 

64x32 

4.64E-09 

-3.86 

8.53E-08 

-3.27 

II 

o 

04x02 

1.02E-05 


7.80E-05 


08x04 

4.72E-07 

-4.43 

4.65E-06 

-4.06 

16x08 

2.14E-08 

-4.46 

2.43E-07 

-4.25 

32x16 

1.15E-09 

-4.22 

1.33E-08 

-4.19 

64x32 

6.83E-11 

-4.06 

7.95E-10 

-4.06 


Table 4 Error convergence is shown for the Navier-Stokes equation. 




Uniform Grid 



Nonuniform Grid 


P=01 

L z error 

L z rate 

L°° error 

L°° rate 

L z error 

L z rate 

L°° error 

L°° rate 

4 

1.95E-01 


5.92E-01 


4.14E-01 


1.19E-00 


8 

7.11E-02 

1.46 

2.73E-01 

1.11 

1.81E-01 

1.19 

7.30E-01 

0.71 

16 

2.10E-02 

1.76 

9.51E-02 

1.53 

5.26E-02 

1.78 

2.52E-01 

1.53 

32 

5.57E-03 

1.92 

2.57E-02 

1.89 

1.51E-02 

1.80 

8.85E-02 

1.51 

64 

1.42E-03 

1.97 

6.50E-03 

1.98 

4.05E-03 

1.91 

2.40E-02 

1.88 

128 

3.58E-04 

1.99 

1.63E-03 

1.99 

1.03E-03 

1.97 

6.06E-03 

1.98 

P=02 

4 

2.64E-02 


7.87E-02 


1.02E-01 


4.11E-01 


8 

6.16E-03 

2.10 

3.04E-02 

1.37 

1.22E-02 

3.07 

5.04E-02 

3.03 

16 

8.22E-04 

2.91 

4.13E-03 

2.88 

5.54E-03 

1.14 

3.74E-02 

0.43 

32 

9.90E-05 

3.05 

8.73E-04 

2.24 

7.46E-04 

2.89 

7.41E-03 

2.33 

64 

1.22E-05 

3.02 

1.09E-04 

2.99 

9.17E-05 

3.03 

8.78E-04 

3.08 

dd 

II 

o 

CO 

4 

8.08E-03 


4.73E-02 


1.54E-02 


3.60E-02 


8 

4.94E-04 

4.03 

3.93E-03 

3.59 

7.22E-03 

1.09 

5.53E-02 

-0.62 

16 

3.45E-05 

3.84 

3.77E-04 

3.38 

3.86E-04 

4.22 

4.03E-03 

3.78 

32 

2.12E-06 

4.03 

2.47E-05 

3.93 

2.89E-05 

3.74 

4.00E-04 

3.34 

64 

1.28E-07 

4.05 

1.46E-06 

4.08 

1.76E-06 

4.04 

2.58E-05 

3.95 

P=04 

4 

1.33E-03 


7.03E-03 


1.16E-02 


8.84E-02 


8 

8.24E-05 

4.02 

8.16E-04 

3.11 

1.12E-03 

3.37 

7.55E-03 

3.55 

16 

2.35E-06 

5.13 

3.21E-05 

4.67 

7.72E-05 

3.86 

8.05E-04 

3.23 

32 

7.15E-08 

5.04 

1.09E-06 

4.88 

2.23E-06 

5.11 

3.13E-05 

4.69 

64 

2.21E-09 

5.02 

3.43E-08 

4.99 

6.92E-08 

5.01 

1.10E-06 

4.83 
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Table A1 Differentiation operators for polynomials of degree one through four. 


p = 

1 

P = 

2 

P = 

3 

P = 

4 

xl = 

-1 

xl = 

-1 

xl = 

-1 

xl = 

-1 

x2 = 

+ 1 

x2 = 

0 

x2 = 

-1/V5 

x2 = 

~\IW 



x3 = 

+1 

x3 = 

+1 / VE 

x3 = 

0 





x4 = 

+i 

x4 = 

+VV7 







x5 = 

+1 

pi = 

+ 1 

pi = 

1/3 

pi = 

1/6 

pi = 

1 / 10 

P 2 = 

+ 1 

P 2 = 

4/3 

P 2 = 

5/6 

P 2 = 

49 / 90 



p3 = 

1/3 

p3 = 

5/6 

p3 = 

32 / 45 





P 4 = 

1/6 

P 4 = 

49 / 90 







p5 = 

1 / 10 

qll = 

-1/2 

qll = 

-1/2 

qll = 

-1/2 

qll = 

-1/2 

ql2 = 

+1/2 

ql2 = 

2/3 

ql2 = 

(5/24 (1 + V5)) 

ql2 = 

(7/120) (7+v^I) 

q22 = 

+1/2 

ql3 = 

-1/6 

ql3 = 

(-5/24 (-1 + x/5)) 

ql3 = 

-(4/15 ) 



q22 = 

0 

ql4 = 

1/12 

ql4 = 

-(7/120) (-7+V21) 



q23 = 

2/3 

q22 = 

0 

ql5 = 

- 1/20 



q33 = 

+1/2 

q23 = 

(5 y/5/12) 

q22 = 

0 





q24 = 

(5/24 (1 - V5)) 

q23 = 

(28 i/773)/ 45 





q33 = 

0 

q24 = 

-(49 i/7/3)/180 





q34 = 

(5/24 (1 + x/5)) 

q25 = 

-(7/120) (-7+V2I) 





q44 = 

1/2 

q33 = 

0 







q34 = 

(28 (V7/3))/45 







q35 = 

- 4/15 







q44 = 

0 







q45 = 

(7/120) (7+v^I) 







q55 = 

1/2 
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